Setup, and Packages

Code Packages:

# Load packages
library(gridExtra) # For using grid.table to save tables as images
## Warning: package 'gridExtra' was built under R version 4.4.3
library(flextable) # Alternative method to save df as image
## Warning: package 'flextable' was built under R version 4.4.3
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'lubridate' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine() masks gridExtra::combine()
## ✖ purrr::compose() masks flextable::compose()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(zoo)
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step
library(rlang)
## 
## Attaching package: 'rlang'
## 
## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice
library(nortest) # For ad.test
library(car) # For homogeneity of variance with leveneTest
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(rcompanion) # For scheirerRayHare test
## Warning: package 'rcompanion' was built under R version 4.4.3
library(pastecs) # For stat.desc
## Warning: package 'pastecs' was built under R version 4.4.3
## 
## Attaching package: 'pastecs'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
library(ggpubr) # For assumption plots
## Warning: package 'ggpubr' was built under R version 4.4.3
## 
## Attaching package: 'ggpubr'
## 
## The following objects are masked from 'package:flextable':
## 
##     border, font, rotate
library(conover.test) # For the conover.Iman() test function
library(DescTools) # Contain alternative ConoverTest function
## Warning: package 'DescTools' was built under R version 4.4.3
## 
## Attaching package: 'DescTools'
## 
## The following object is masked from 'package:car':
## 
##     Recode
library(rstatix) # tidyverse adapted stat tests and mshapiro_test() functions
## Warning: package 'rstatix' was built under R version 4.4.3
## 
## Attaching package: 'rstatix'
## 
## The following object is masked from 'package:stats':
## 
##     filter
library(DHARMa) # For multi-level lm models and shapiro test
## Warning: package 'DHARMa' was built under R version 4.4.3
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(patchwork) # For combining ggplots
## Warning: package 'patchwork' was built under R version 4.4.3

Data Entry: Read-in data

# Read and call data into df

# Primary samples df
ABL90 <- read.csv("Raw_Data/ABL90_Raw.csv")

# River's Parturition metadata
partuition_subcat <- read.csv("Raw_Data/River's_Rockfish_Metadata_Parturition_V7.csv")

Data Filtering

# Data sifting: ABL90 dataset

# Step 1: Remove missing info
ABL_set1 <- ABL90 %>%
  filter(Patient.ID_edited != "") 
  

# Step 2: Separate Patient.ID into columns of sample types: blood plasma and instant freeze plasma 
ABL_set2 <- separate(ABL_set1, 'Patient.ID_edited', into = c("Patient.ID_edited", "Sample_Type"), sep = ",")  
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 4 rows [489, 619, 1034,
## 1177].
# Step 3:
# Convert Patient.ID_edited column data to character type
ABL_set2$Patient.ID_edited <- as.character(ABL_set2$Patient.ID_edited) 
# Step 4 & 5: Sussing out specific sample errors  

# Remove insufficient samples
ABL_set3 <- ABL_set2 %>%
  filter(!is.na("Type")) %>%
  filter(!str_detect(Errors.detected.during.measurement, "Insufficient sample"))

# Remove inhomogeneous samples
ABL_set3 <- ABL_set3 %>%
  filter(!is.na("Type")) %>%
  filter(!str_detect(Errors.detected.during.measurement, "Inhomogeneous sample"))

# Step 6: Filter for only blood samples
ABL_b_samp <- ABL_set3 %>%
  filter(Sample_Type == "b")

Merge Datasets:

Join ABL90 df with parturition_subcat df: Making ABL_merged

# Rename columns to match: Change 'ID' in parturition_subcat to 'Patient.ID_edited' so it is ready to join with ABL90 df

# Please note the 'Treatment' col in parturition_subcat does not match with the finalized 'Ambient.Or.OAH' col in metadata_atresia_guide or ABL90_merged df

# Rename 'ID' col to 'patient.ID_edited'
partuition_subcat <- partuition_subcat %>%
  rename(Patient.ID_edited = ID)
# Connect main ABL90 dataset with my (River's) sub categorized parturition metadata 
# ABL_merged <- partuition_subcat %>%
# left_join(ABL_b_samp, by = "Patient.ID_edited")

ABL_merged <- ABL_b_samp %>%
inner_join(partuition_subcat, by = "Patient.ID_edited")

Rename Columns Neatly:

# Rename Treatment & Parturition Outcome
ABL_merged <- ABL_merged %>%
  rename(Ambient.Or.OAH = Consensus_General_Treatment,
         Pregnant.Or.Atresia = Consensus_Atresia_Or_Pregnant)

Remove Samples: Mortalities, NA’s, or Missing Info, and Replicates

Remove Replicate ID’s

# Cross Check all Columns: Cross validation to suss correct replicate row
check_params <- ABL_merged %>%
  select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
  filter(Patient.ID_edited == "9782D")
print(check_params)
##   Patient.ID_edited             Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1             9782D 12/28/2023 10:04     1016         C 65uL        Ambient
## 2             9782D   1/9/2024 17:21      863         S 65uL        Ambient
##   Pregnant.Or.Atresia    pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1    Post Parturition 7.348          0.5     9.8          184          149
## 2    Post Parturition 7.446          0.7    13.3          181          161
##   K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1         2.8          1.32         4.2
## 2         2.9          1.31         6.8
check_params <- ABL_merged %>%
  select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
  filter(Patient.ID_edited == "9783D")
print(check_params)
##   Patient.ID_edited            Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1             9783D 1/19/2024 10:37      934         S 65uL        Ambient
## 2             9783D  1/7/2024 12:32      860         S 65uL        Ambient
##   Pregnant.Or.Atresia    pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1    Post Parturition 7.484          1.7     0.2          175          193
## 2    Post Parturition 7.450          1.3    14.7          179          164
##   K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1        34.5          1.29         8.5
## 2         2.5          1.30         5.3
# Removed replicates: 9782D (2x) and 9783D (2x)
ABL_merged <- ABL_merged %>%
  filter(Sample.. != "863",
         Sample.. != "934")

Review Replicates:

Row ID Time Sample.. Measuring.Mode pH Glu.mmol.L. Hct…. Na…mmol.L. Cl…mmol.L. K…mmol.L. Ca.mmol.L.
32 9782D 12/28/2023 10:04 1016 C 65uL 7.348 0.5 9.8 184 149 2.8 1.32
39 9782D 1/9/2024 17:21 863 S 65uL 7.446 0.7 13.3 181 161 2.9 1.31
36 9783D 1/19/2024 10:37 934 S 65uL 7.484 1.7 0.2 175 193 34.5 1.29
41 9783D 1/7/2024 12:32 860 S 65uL 7.450 1.3 14.7 179 164 2.5 1.30

Methods for replicate sample removal: Check processing date of IDs, and remove samples that do not match that date.

Case ID_9782D:

Deciding info: - For ID-9782D, refer to metadata sheet and look for other samples with that same processing date, then view thoes samples ABL90 ‘Time’. Try to find a matching ‘Time’ with other IDs and 9782D that were processed together during the same time recorded by the machine.

Case ID_9783D

Samples: Data subsets

Remove Motalities and No info IDs and assign analysis ready data-frames:

# New df with Moralities removed: Note none of these samples made it into ABL90 df anyway, so looks like they are already filtered out.
 Live_Samples <- ABL_merged %>%
  filter(Patient.ID_edited != "9780C", # Mortality
         Patient.ID_edited != "777AE", # Mortality
         Patient.ID_edited != "777CA") # Mortality after parturition


# New df with mortality and 'No info' Id's removed
 Primary_Samples <- Live_Samples %>%
   filter(Patient.ID_edited != "777A0", # No info
           Patient.ID_edited != "9782F", # No info
           Patient.ID_edited != "777B3", # No info
           Patient.ID_edited != "777AA", # No info
           Patient.ID_edited != "777DE", # No info
           Patient.ID_edited != "777CE", # No info
           Patient.ID_edited != "777A6") # No info
 
 
# New df of Only Ambient Treatment: For testing parturition success
 Ambient_Samples <- Primary_Samples %>%
   filter(Ambient.Or.OAH == "Ambient")
 
# New df of Fecundity Samples:
 Fecundity_Samples <- Primary_Samples %>%
   filter(Fecundity_Or_Brood_Count != "NA",
     Fecundity_Class != "NA") # removes 97706

# Fecundity samples without atresia Ids
 Fecundity_No_Atresia_Samples <- Primary_Samples %>%
   filter(All_Fecundity != "Na",
          All_Fecundity != 0)  # removes atresia IDs

# Fecundity Ambient Samples
 Amb_Fec_Samples <- Primary_Samples %>%
   filter(Ambient.Or.OAH == "Ambient",
          All_Fecundity != "NA") # removes OAH treatment samples and ID 97706
 
# Fecundity Ambient Samples without atresia Ids
 Amb_Fec_No_Atresia_Samples <- Primary_Samples %>%
   filter(Ambient.Or.OAH == "Ambient", # removes OAH treatment 
          All_Fecundity != "NA", # removes all missing info IDs
          All_Fecundity != 0) # removes atresia IDs

Change Data Types:

# Change Continuous Data to type to numeric, double, or factor

# Change data type of ions to factor or numeric
#glimpse(General_Samples)
#glimpse(Ambient_Samples)
 
Primary_Samples <- Primary_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Ambient_Samples <- Ambient_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Fecundity_Samples <- Fecundity_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Fecundity_No_Atresia_Samples <- Fecundity_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Amb_Fec_Samples <- Amb_Fec_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
  
Amb_Fec_No_Atresia_Samples <- Amb_Fec_No_Atresia_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

# glimpse(Primary_Samples)
# glimpse(Ambient_Samples)
# glimpse(Fecundity_Samples)
# glimpse(Amb_Fec_Samples)

Order Factor Levels

# Fecundity data 

#unique(Fecundity_Samples$Fecundity_Class)

# Arrange the order of parturition categories for visualizations & tests

# Primary Samples

# Arrange Treatment
Primary_Samples$Ambient.Or.OAH <- ordered(Primary_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))

# Arrange Parturition outcome
Primary_Samples$Pregnant.Or.Atresia <- ordered(Primary_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))

# Arrange Brood Condition
Primary_Samples$Consensus_Brood_Condition <- ordered(Primary_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Very Low", "Atresia"))

# Arrange Fecundity Class
Primary_Samples$Fecundity_Class <- ordered(Primary_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))

# Remove NAs from ordered character factors
Primary_Samples <- Primary_Samples %>%
  filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
  filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))
 
 
# Ambient data
Ambient_Samples$Ambient.Or.OAH <- ordered(Ambient_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))

Ambient_Samples$Pregnant.Or.Atresia <- ordered(Ambient_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))

Ambient_Samples$Consensus_Brood_Condition <- ordered(Ambient_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Atresia"))

Ambient_Samples$Fecundity_Class <- ordered(Ambient_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))

# Remove NAs from ordered character factors
Ambient_Samples <- Ambient_Samples %>%
  filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
  filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))

# Ambient Fecundity data
Amb_Fec_No_Atresia_Samples$Ambient.Or.OAH <- ordered(Amb_Fec_No_Atresia_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))

Amb_Fec_No_Atresia_Samples$Pregnant.Or.Atresia <- ordered(Amb_Fec_No_Atresia_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))

Amb_Fec_No_Atresia_Samples$Consensus_Brood_Condition <- ordered(Amb_Fec_No_Atresia_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Atresia"))

Amb_Fec_No_Atresia_Samples$Fecundity_Class <- ordered(Amb_Fec_No_Atresia_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))

# Remove NAs from ordered character factors
Amb_Fec_No_Atresia_Samples <- Amb_Fec_No_Atresia_Samples %>%
  filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
  filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))

Save filtered dataset into data-worked folder

# Save merged dataset in a worked folder

write.csv(x = Primary_Samples, file = "Worked-Data/Primary_Samples")

write.csv(x = Ambient_Samples, file = "Worked-Data/Ambient_Samples")

write.csv(x = Fecundity_Samples, file = "Worked-Data/Fecundity_Samples")

write.csv(x = Fecundity_No_Atresia_Samples, file = "Worked-Data/Fecundity_No_Atresia_Samples")

write.csv(x = Amb_Fec_Samples, file = "Worked-Data/Amb-Fec_Samples")

write.csv(x = Amb_Fec_No_Atresia_Samples, file = "Worked-Data/Amb_Fec_No_Atresia_Samples")

Sample Size: n

# Sample Size

# Ambient data:
ambient_sample_size_table <- Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(ambient_sample_size_table)
## # A tibble: 4 × 3
##   Ambient.Or.OAH Consensus_Brood_Condition n_size
##   <ord>          <ord>                      <int>
## 1 Ambient        Excellent                      1
## 2 Ambient        Normal                         4
## 3 Ambient        Low                            3
## 4 Ambient        Atresia                       13
 pdf("data-images/ambient_sample_size_table.pdf")
 grid.table(ambient_sample_size_table)
 dev.off()
## png 
##   2
 png("data-images/ambient_sample_size_table.png")
 grid.table(ambient_sample_size_table)
 dev.off()
## png 
##   2
ambient_fecundity_class_sample_size_table <- Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(ambient_fecundity_class_sample_size_table)
## # A tibble: 3 × 3
##   Ambient.Or.OAH Fecundity_Class n_size
##   <ord>          <ord>            <int>
## 1 Ambient        High (>50,000)       4
## 2 Ambient        Low (~1,000s)        4
## 3 Ambient        Atresia             13
 pdf("data-images/ambient_fecundity_class_sample_size_table.pdf")
 grid.table(ambient_fecundity_class_sample_size_table)
 dev.off()
## png 
##   2
 png("data-images/ambient_fecundity_class_sample_size_table.png")
 grid.table(ambient_fecundity_class_sample_size_table)
 dev.off()
## png 
##   2

Sample Size Barplot:

# View Sample Distributions
# Sample Size barplot
# Ambient Samples

# Consensus brood condition sample barplot
ambient_sample_size_barplot <- Ambient_Samples %>%
  group_by(Consensus_Brood_Condition, Ambient.Or.OAH) %>% 
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
  ggplot(aes(x = Consensus_Brood_Condition, y = n_size)) +
  geom_bar(stat = "identity") +
  facet_grid(~ Ambient.Or.OAH) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
  labs(title = "Ambient Sample Size",
        x = "Parturition Outcome",
        y = "Sample Size (n)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(ambient_sample_size_barplot)

ggsave(filename = "data-images/.png", plot = ambient_sample_size_barplot, device = "png")
## Saving 7 x 5 in image
ggsave(filename = "data-images/ambient_sample_size_barplot.pdf", plot = ambient_sample_size_barplot, device = "pdf")
## Saving 7 x 5 in image
# Fecundity class sample barplot
ambient_sample_size_barplot2 <- Ambient_Samples %>%
  group_by(Fecundity_Class, Ambient.Or.OAH) %>% 
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
  ggplot(aes(x = Fecundity_Class, y = n_size)) +
  geom_bar(stat = "identity") +
  facet_grid(~ Ambient.Or.OAH) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
  labs(title = "Ambient Sample Size",
        x = "Parturition Outcome",
        y = "Sample Size (n)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(ambient_sample_size_barplot2)

ggsave(filename = "data-images/.png", plot = ambient_sample_size_barplot2, device = "png")
## Saving 7 x 5 in image
ggsave(filename = "data-images/ambient_sample_size_barplot2.pdf", plot = ambient_sample_size_barplot2, device = "pdf")
## Saving 7 x 5 in image
# Amb Fec No atresia sample size

amb_fec_no_atresia_barplot <- Amb_Fec_No_Atresia_Samples %>%
  group_by(Fecundity_Class, Ambient.Or.OAH) %>% 
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
  ggplot(aes(x = Fecundity_Class, y = n_size)) +
  geom_bar(stat = "identity") +
  facet_grid(~ Ambient.Or.OAH) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(4, 8, 12)) +
  labs(title = "Ambient No Atresia Sample Size",
        x = "Parturition Outcome",
        y = "Sample Size (n)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(amb_fec_no_atresia_barplot)

ggsave(filename = "data-images/.png", plot = amb_fec_no_atresia_barplot, device = "png")
## Saving 7 x 5 in image
ggsave(filename = "data-images/amb_fec_no_atresia_barplot.pdf", plot = amb_fec_no_atresia_barplot, device = "pdf")
## Saving 7 x 5 in image

Paramater Analysis:

pH

pH Summary Stats

# pH
# Summary stats

# Ambient data
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
  summarize(count = n(),
            median = round(median(pH), 3),
            mean = round(mean(pH), 3),
            sd = round(sd(pH), 3),
            cv = round(sd(pH)/mean(pH), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
##   Ambient.Or.OAH Consensus_Brood_Condition count median  mean     sd     cv
##   <ord>          <ord>                     <int>  <dbl> <dbl>  <dbl>  <dbl>
## 1 Ambient        Excellent                     1   7.43  7.43 NA     NA    
## 2 Ambient        Normal                        4   7.48  7.45  0.087  0.012
## 3 Ambient        Low                           3   7.43  7.44  0.046  0.006
## 4 Ambient        Atresia                      13   7.37  7.38  0.1    0.014
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(pH), 3),
            mean = round(mean(pH), 3),
            sd = round(sd(pH), 3),
            cv = round(sd(pH)/mean(pH), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   7.44  7.43 0.076 0.01 
## 2 Ambient        Low (~1,000s)       4   7.46  7.45 0.054 0.007
## 3 Ambient        Atresia            13   7.37  7.38 0.1   0.014
Amb_Fec_No_Atresia_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(pH), 3),
            mean = round(mean(pH), 3),
            sd = round(sd(pH), 3),
            cv = round(sd(pH)/mean(pH), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   7.44  7.43 0.076 0.01 
## 2 Ambient        Low (~1,000s)       4   7.46  7.45 0.054 0.007

pH Boxplot

# pH boxplot

# Ambient Samples
pH_ambient_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = pH)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = pH), alpha = 0.5, colour = "black") +
  labs(title = "pH", x = "Parturition Success", y = "pH") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_ambient_boxplot)

ggsave(filename = "data-images/pH_ambient_boxplot.pdf", plot = pH_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/pH_ambient_boxplot.png", plot = pH_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
pH_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = pH)) +
  geom_point(aes(x = Fecundity_Class, y = pH), alpha = 0.5, colour = "black") +
  labs(title = "pH", x = "Fecundity Class", y = "pH") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_ambient_fecundity_boxplot)

ggsave(filename = "data-images/pH_ambient_fecundity_boxplot.pdf", plot = pH_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/pH_ambient_fecundity_boxplot.png", plot = pH_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
pH_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = pH)) +
  geom_point(aes(x = Fecundity_Class, y = pH), alpha = 0.5, colour = "black") +
  labs(title = "pH", x = "Fecundity Class", y = "pH") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_amb_fec_no_atresia_boxplot)

ggsave(filename = "data-images/pH_amb_fec_no_atresia_boxplot.pdf", plot = pH_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/pH_amb_fec_no_atresia_boxplot.png", plot = pH_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
pH_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH), colour = "black") +
  labs(title = "pH", x = "Weight Adjusted Fecundity", y = "pH") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_amb_fec_no_atresia_scatterplot)

ggsave(filename = "data-images/pH_amb_fec_no_atresia_scatterplot.pdf", plot = pH_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/pH_amb_fec_no_atresia_scatterplot.png", plot = pH_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image

pH Stat Tests

Differences: Between Brood Conditon

# pH
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

pH_aov_ambient_brood <- aov(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(pH_aov_ambient_brood)
##                           Df  Sum Sq  Mean Sq F value Pr(>F)
## Consensus_Brood_Condition  3 0.02079 0.006930   0.795  0.513
## Residuals                 17 0.14815 0.008715
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(pH_aov_ambient_brood)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
## $Consensus_Brood_Condition
##                          diff        lwr        upr     p adj
## Normal-Excellent   0.01675000 -0.2799312 0.31343123 0.9984685
## Low-Excellent      0.00400000 -0.3024111 0.31041106 0.9999809
## Atresia-Excellent -0.05430769 -0.3296845 0.22106915 0.9423449
## Low-Normal        -0.01275000 -0.2154219 0.18992187 0.9978869
## Atresia-Normal    -0.07105769 -0.2227829 0.08066756 0.5566650
## Atresia-Low       -0.05830769 -0.2282740 0.11165858 0.7651322
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$pH)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$pH
## W = 0.98534, p-value = 0.9808
shapiro.test(residuals(aov(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(pH ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.95636, p-value = 0.446
# QQplot:
pH_ambient_res_qqplot <- ggqqplot(residuals(aov(pH ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Ambient pH Residual QQplot",
     subtitle = "residuals(aov(pH ~ Consensus_Brood_Condition, data = Ambient_Samples))",
                               x = "pH Theoretical Quantiles (Predicted values)",
                               y = "pH Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(pH_ambient_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
#bartlett.test(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)

# Nonparametric variance test (Levene's test):
leveneTest(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   0.671 0.5815
##       17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  pH by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 2.8578, df = 3, p-value = 0.4141
# Nonparametric Post Hoc: Dunn's Test
DunnettTest(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $Excellent
##                          diff     lwr.ci    upr.ci   pval    
## Normal-Excellent   0.01675000 -0.2376275 0.2711275 0.9925    
## Low-Excellent      0.00400000 -0.2587200 0.2667200 0.9999    
## Atresia-Excellent -0.05430769 -0.2904186 0.1818032 0.8155    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pH Ambient Results

Which brood condition significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): - Variance (Bartletts test): NA

Differences: Between Fecundity Class

# pH
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

pH_aov_ambient_fec <- aov(pH ~ Fecundity_Class, data = Ambient_Samples)
summary(pH_aov_ambient_fec)
##                 Df  Sum Sq  Mean Sq F value Pr(>F)
## Fecundity_Class  2 0.02183 0.010916   1.336  0.288
## Residuals       18 0.14711 0.008173
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(pH_aov_ambient_fec)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = pH ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                     diff        lwr        upr     p adj
## Low (~1,000s)-High (>50,000)  0.02675000 -0.1363959 0.18989587 0.9084726
## Atresia-High (>50,000)       -0.05080769 -0.1827287 0.08111329 0.5966404
## Atresia-Low (~1,000s)        -0.07755769 -0.2094787 0.05436329 0.3141611
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$pH)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$pH
## W = 0.98534, p-value = 0.9808
shapiro.test(residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.96555, p-value = 0.6339
# QQplot:
pH_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient pH Residual QQplot",
     subtitle = "residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "pH Theoretical Quantiles (Predicted values)",
                               y = "pH Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(pH_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(pH ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  pH by Fecundity_Class
## Bartlett's K-squared = 1.4521, df = 2, p-value = 0.4838
# Nonparametric variance test (Levene's test):
leveneTest(pH ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2    0.44 0.6508
##       18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(pH ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  pH by Fecundity_Class
## Kruskal-Wallis chi-squared = 3.0375, df = 2, p-value = 0.219
# Nonparametric Post Hoc: Dunn's test
DunnettTest(pH ~ Fecundity_Class, data = Ambient_Samples)
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`High (>50,000)`
##                                     diff     lwr.ci     upr.ci   pval    
## Low (~1,000s)-High (>50,000)  0.02675000 -0.1252448 0.17874478 0.8717    
## Atresia-High (>50,000)       -0.05080769 -0.1737118 0.07209643 0.5071    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pH Ambient Results

Which fecundity class significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): Met

Nonparametric variance test (Levene’s test):

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# Ph
# Linear Model
pH_amb_fec_no_atresia_lm <- lm(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(pH_amb_fec_no_atresia_lm)
## 
## Call:
## lm(formula = pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Fecundity_No_Atresia_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.11940 -0.02163  0.01153  0.04429  0.05591 
## 
## Coefficients:
##                                                   Estimate Std. Error t value
## (Intercept)                                      7.4752402  0.0369759  202.16
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -0.0001977  0.0001939   -1.02
##                                                 Pr(>|t|)    
## (Intercept)                                     4.01e-16 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.338    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05852 on 8 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.115,  Adjusted R-squared:  0.004385 
## F-statistic:  1.04 on 1 and 8 DF,  p-value: 0.3377
pH_ambient_sim <- simulateResiduals(fittedModel = pH_amb_fec_no_atresia_lm) 
plot(pH_ambient_sim)
## qu = 0.25, log(sigma) = -4.135132 : outer Newton did not converge fully.

plotResiduals(pH_amb_fec_no_atresia_lm)
## qu = 0.25, log(sigma) = -4.135132 : outer Newton did not converge fully.

testDispersion(pH_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Normality test
shapiro.test(residuals(pH_amb_fec_no_atresia_lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(pH_amb_fec_no_atresia_lm)
## W = 0.89785, p-value = 0.2075
# Transform data: square root (undo log scale?)
pH_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(pH = sqrt(pH))

# Rerun lm model with square root transformed data
pH_amb_fec_no_atresia_lm_sqrt <- lm(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = pH_ambient_sqrt)
summary(pH_amb_fec_no_atresia_lm_sqrt)
## 
## Call:
## lm(formula = pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = pH_ambient_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.020889 -0.004920  0.002934  0.007526  0.010787 
## 
## Coefficients:
##                                                   Estimate Std. Error t value
## (Intercept)                                      2.734e+00  7.412e-03 368.821
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -4.078e-05  4.222e-05  -0.966
##                                                 Pr(>|t|)    
## (Intercept)                                     2.68e-14 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.371    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01157 on 6 degrees of freedom
## Multiple R-squared:  0.1346, Adjusted R-squared:  -0.009642 
## F-statistic: 0.9331 on 1 and 6 DF,  p-value: 0.3714
pH_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = pH_amb_fec_no_atresia_lm_sqrt) 
plot(pH_amb_fec_sqrt_sim)
## qu = 0.25, log(sigma) = -2.642951 : outer Newton did not converge fully.

plotResiduals(pH_amb_fec_no_atresia_lm_sqrt)
## qu = 0.25, log(sigma) = -2.642951 : outer Newton did not converge fully.

testDispersion(pH_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
pH_amb_fec_no_atresia_lm_res <- residuals(pH_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(pH_amb_fec_no_atresia_lm_res) # normal
## 
##  Shapiro-Wilk normality test
## 
## data:  pH_amb_fec_no_atresia_lm_res
## W = 0.90414, p-value = 0.3147
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)

Residuals Plot (if significant)

# pH
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

pH_ambient_lm_scatterplot <- ggplot(data = pH_ambient_sqrt, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Weight Adjusted Fecudnity", y = "Sqrt pH", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/pH_ambient_lm_scatterplot.pdf", plot = pH_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/pH_ambient_lm_scatterplot.png", plot = pH_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# sodium
# lm Boxplot

pH_ambient_lm_boxplot <- ggplot(data = pH_ambient_sqrt, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Weight Adjusted Fecundity", y = "Sqrt pH", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(pH_ambient_lm_boxplot)
## Warning: Continuous x aesthetic
## ℹ did you forget `aes(group = ...)`?
## `geom_smooth()` using formula = 'y ~ x'

Hct : Hematocrit

hct Summary Stats

# hct
# Summary stats

# Ambient data
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
  summarize(count = n(),
            median = round(median(Hct....), 3),
            mean = round(mean(Hct....), 3),
            sd = round(sd(Hct....), 3),
            cv = round(sd(Hct....)/mean(Hct....), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
##   Ambient.Or.OAH Consensus_Brood_Condition count median  mean    sd     cv
##   <ord>          <ord>                     <int>  <dbl> <dbl> <dbl>  <dbl>
## 1 Ambient        Excellent                     1   18.2  18.2 NA    NA    
## 2 Ambient        Normal                        4   13.8  15.4  4.33  0.281
## 3 Ambient        Low                           3   14    13.4  3.29  0.245
## 4 Ambient        Atresia                      13   17.3  18.2  2.33  0.128
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(Hct....), 3),
            mean = round(mean(Hct....), 3),
            sd = round(sd(Hct....), 3),
            cv = round(sd(Hct....)/mean(Hct....), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   16.4  16.7  4.12 0.247
## 2 Ambient        Low (~1,000s)       4   13.4  13.3  2.70 0.203
## 3 Ambient        Atresia            13   17.3  18.2  2.33 0.128
Amb_Fec_No_Atresia_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(Hct....), 3),
            mean = round(mean(Hct....), 3),
            sd = round(sd(Hct....), 3),
            cv = round(sd(Hct....)/mean(Hct....), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   16.4  16.7  4.12 0.247
## 2 Ambient        Low (~1,000s)       4   13.4  13.3  2.70 0.203

hct Boxplot

# hct boxplot

# Ambient Samples
hct_ambient_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = Hct....)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = Hct....), alpha = 0.5, colour = "black") +
  labs(title = "Hematocrit", x = "Parturition Success", y = "Hematocrit (%)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(hct_ambient_boxplot)

ggsave(filename = "data-images/hct_ambient_boxplot.pdf", plot = hct_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/hct_ambient_boxplot.png", plot = hct_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
hct_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Hct....)) +
  geom_point(aes(x = Fecundity_Class, y = Hct....), alpha = 0.5, colour = "black") +
  labs(title = "Hematocrit", x = "Fecundity Class", y = "Hematocrit (%)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(hct_ambient_fecundity_boxplot)

ggsave(filename = "data-images/hct_ambient_fecundity_boxplot.pdf", plot = hct_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/hct_ambient_fecundity_boxplot.png", plot = hct_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
hct_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Hct....)) +
  geom_point(aes(x = Fecundity_Class, y = Hct....), alpha = 0.5, colour = "black") +
  labs(title = "Hematocrit", x = "Fecundity Class", y = "Hematocrit (%)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(hct_amb_fec_no_atresia_boxplot)

ggsave(filename = "data-images/hct_amb_fec_no_atresia_boxplot.pdf", plot = hct_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/hct_amb_fec_no_atresia_boxplot.png", plot = hct_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
hct_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....), colour = "black") +
  labs(title = "Hematocrit", x = "Weight Adjusted Fecundity", y = "Hematocrit (%)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(hct_amb_fec_no_atresia_scatterplot)

ggsave(filename = "data-images/hct_amb_fec_no_atresia_scatterplot.pdf", plot = hct_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/hct_amb_fec_no_atresia_scatterplot.png", plot = hct_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image

Hct Stat Tests

Differences: Between Brood Conditon

# hct
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

hct_aov_ambient_brood <- aov(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(hct_aov_ambient_brood) # Significant
##                           Df Sum Sq Mean Sq F value Pr(>F)  
## Consensus_Brood_Condition  3  68.88  22.959   2.735 0.0758 .
## Residuals                 17 142.69   8.394                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(hct_aov_ambient_brood) # Atresia-Low Significant
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
## $Consensus_Brood_Condition
##                            diff         lwr       upr     p adj
## Normal-Excellent  -2.825000e+00 -12.0325283  6.382528 0.8190431
## Low-Excellent     -4.766667e+00 -14.2761610  4.742828 0.5017753
## Atresia-Excellent  3.552714e-15  -8.5463446  8.546345 1.0000000
## Low-Normal        -1.941667e+00  -8.2316060  4.348273 0.8163755
## Atresia-Normal     2.825000e+00  -1.8838065  7.533807 0.3513279
## Atresia-Low        4.766667e+00  -0.5082517 10.041585 0.0843920
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Hct....)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Hct....
## W = 0.97482, p-value = 0.8355
shapiro.test(residuals(aov(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.9585, p-value = 0.4864
# QQplot:
hct_ambient_res_qqplot <- ggqqplot(residuals(aov(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Ambint Hematocrit Residual QQplot",
     subtitle = "residuals(aov(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples))",
                               x = "Hct Theoretical Quantiles (Predicted values)",
                               y = "Hct Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(hct_ambient_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)

# Nonparametric variance test (Levene's test):
leveneTest(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   0.507 0.6826
##       17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Hct.... by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 6.6608, df = 3, p-value = 0.08353
# Nonparametric Post Hoc: Dunn's test
DunnettTest(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $Excellent
##                        diff     lwr.ci   upr.ci   pval    
## Normal-Excellent  -2.825000 -10.719629 5.069629 0.5993    
## Low-Excellent     -4.766667 -12.920205 3.386871 0.2885    
## Atresia-Excellent  0.000000  -7.327724 7.327724 1.0000    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hematocrit Ambient Results

Which brood condition significantly different (if any)? - Atresia-Low (Tukey HSD: p adj = 0.0843920)

Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartletts test): NA

Differences: Between Fecundity Class

# hct
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

hct_aov_ambient_fec <- aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
summary(hct_aov_ambient_fec)
##                 Df Sum Sq Mean Sq F value Pr(>F)  
## Fecundity_Class  2  73.83   36.92   4.824  0.021 *
## Residuals       18 137.74    7.65                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(hct_aov_ambient_fec)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Hct.... ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                              diff        lwr      upr     p adj
## Low (~1,000s)-High (>50,000) -3.4 -8.3921454 1.592145 0.2186627
## Atresia-High (>50,000)        1.5 -2.5366864 5.536686 0.6176987
## Atresia-Low (~1,000s)         4.9  0.8633136 8.936686 0.0162654
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Hct....)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Hct....
## W = 0.97482, p-value = 0.8355
shapiro.test(residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.97094, p-value = 0.7538
# QQplot:
hct_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Hematocrit Residual QQplot",
     subtitle = "residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Hct Theoretical Quantiles (Predicted values)",
                               y = "Hct Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(hct_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Hct.... by Fecundity_Class
## Bartlett's K-squared = 1.7219, df = 2, p-value = 0.4228
# Nonparametric variance test (Levene's test):
leveneTest(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.2003 0.3241
##       18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Hct.... by Fecundity_Class
## Kruskal-Wallis chi-squared = 6.7821, df = 2, p-value = 0.03367
# Nonparametric Post Hoc: Dunn's Test
DunnettTest(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`High (>50,000)`
##                              diff    lwr.ci   upr.ci   pval    
## Low (~1,000s)-High (>50,000) -3.4 -8.050930 1.250930 0.1642    
## Atresia-High (>50,000)        1.5 -2.260777 5.260777 0.5289    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hct Ambient Results

Which fecundity class significantly different (if any)? - Atresia-Low (~1,000s): (Tukey HSD: p adj = 0.0162654)

Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): Met

Nonparametric variance test (Levene’s test): Met

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# hct
# Linear Model
hct_amb_fec_no_atresia_lm <- lm(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(hct_amb_fec_no_atresia_lm)
## 
## Call:
## lm(formula = Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Fecundity_No_Atresia_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4692 -1.9674 -0.8183  1.4654  7.1131 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                     13.09179    2.12136   6.171
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight  0.01313    0.01112   1.181
##                                                 Pr(>|t|)    
## (Intercept)                                     0.000268 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.271681    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.357 on 8 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.1484, Adjusted R-squared:  0.04192 
## F-statistic: 1.394 on 1 and 8 DF,  p-value: 0.2717
hct_ambient_sim <- simulateResiduals(fittedModel = hct_amb_fec_no_atresia_lm) 
plot(hct_ambient_sim)
## qu = 0.5, log(sigma) = -3.130892 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.137393 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.148809 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.14559 : outer Newton did not converge fully.

plotResiduals(hct_amb_fec_no_atresia_lm)
## qu = 0.5, log(sigma) = -3.130892 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.137393 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.148809 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.14559 : outer Newton did not converge fully.

testDispersion(hct_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Normality test
shapiro.test(residuals(hct_amb_fec_no_atresia_lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(hct_amb_fec_no_atresia_lm)
## W = 0.89945, p-value = 0.2161
# Transform data: square root (undo log scale?)
hct_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(Hct.... = sqrt(Hct....))

# Rerun lm model with square root transformed data
hct_amb_fec_no_atresia_lm_sqrt <- lm(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt)
summary(hct_amb_fec_no_atresia_lm_sqrt)
## 
## Call:
## lm(formula = Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = hct_ambient_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45751 -0.18537 -0.07738  0.06130  0.87682 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                     3.549914   0.292111  12.153
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.002034   0.001664   1.223
##                                                 Pr(>|t|)    
## (Intercept)                                     1.89e-05 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.267    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.456 on 6 degrees of freedom
## Multiple R-squared:  0.1995, Adjusted R-squared:  0.06603 
## F-statistic: 1.495 on 1 and 6 DF,  p-value: 0.2673
hct_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = hct_amb_fec_no_atresia_lm_sqrt) 
plot(hct_amb_fec_sqrt_sim)

plotResiduals(hct_amb_fec_no_atresia_lm_sqrt)

testDispersion(hct_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
hct_amb_fec_no_atresia_lm_res <- residuals(hct_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(hct_amb_fec_no_atresia_lm_res) 
## 
##  Shapiro-Wilk normality test
## 
## data:  hct_amb_fec_no_atresia_lm_res
## W = 0.87306, p-value = 0.1615
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt)

Hct Ambient No Atresia Results

Does a significant correlation exist between weight adj fecundity and parameter?

DHARMA Warning: Data dispersion error

Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): NA

Nonparametric variance test (Levene’s test): NA

Residuals Plot (if significant)

# hct
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

hct_ambient_lm_scatterplot <- ggplot(data = hct_ambient_sqrt, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Hematocrit", x = "Weight Adjusted Fecudnity", y = "Sqrt Hematocrit (%)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(hct_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/hct_ambient_lm_scatterplot.pdf", plot = hct_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/hct_ambient_lm_scatterplot.png", plot = hct_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

Glucose

glucose Summary Stats

# glucose
# Summary stats

# Ambient data
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
  summarize(count = n(),
            median = round(median(Glu..mmol.L.), 3),
            mean = round(mean(Glu..mmol.L.), 3),
            sd = round(sd(Glu..mmol.L.), 3),
            cv = round(sd(Glu..mmol.L.)/mean(Glu..mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
##   Ambient.Or.OAH Consensus_Brood_Condition count median  mean     sd     cv
##   <ord>          <ord>                     <int>  <dbl> <dbl>  <dbl>  <dbl>
## 1 Ambient        Excellent                     1   1.8   1.8  NA     NA    
## 2 Ambient        Normal                        4   1.55  1.58  0.275  0.175
## 3 Ambient        Low                           3   1.1   1.6   1.04   0.653
## 4 Ambient        Atresia                      13   2     2.1   0.698  0.332
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(Glu..mmol.L.), 3),
            mean = round(mean(Glu..mmol.L.), 3),
            sd = round(sd(Glu..mmol.L.), 3),
            cv = round(sd(Glu..mmol.L.)/mean(Glu..mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   1.55  1.55 0.238 0.154
## 2 Ambient        Low (~1,000s)       4   1.5   1.68 0.866 0.517
## 3 Ambient        Atresia            13   2     2.1  0.698 0.332
Amb_Fec_No_Atresia_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(Glu..mmol.L.), 3),
            mean = round(mean(Glu..mmol.L.), 3),
            sd = round(sd(Glu..mmol.L.), 3),
            cv = round(sd(Glu..mmol.L.)/mean(Glu..mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   1.55  1.55 0.238 0.154
## 2 Ambient        Low (~1,000s)       4   1.5   1.68 0.866 0.517

Glucose Boxplot

# glucose boxplot

# Ambient Samples
glucose_ambient_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = Glu..mmol.L.)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Glucose", x = "Parturition Success", y = "Glucose (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(glucose_ambient_boxplot)

ggsave(filename = "data-images/glucose_ambient_boxplot.pdf", plot = glucose_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/glucose_ambient_boxplot.png", plot = glucose_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
glucose_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Glu..mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Glucose", x = "Fecundity Class", y = "Glucose (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(glucose_ambient_fecundity_boxplot)

ggsave(filename = "data-images/glucose_ambient_fecundity_boxplot.pdf", plot = glucose_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/glucose_ambient_fecundity_boxplot.png", plot = glucose_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
glucose_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Glu..mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Glucose", x = "Fecundity Class", y = "Glucose (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(glucose_amb_fec_no_atresia_boxplot)

ggsave(filename = "data-images/glucose_amb_fec_no_atresia_boxplot.pdf", plot = glucose_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/glucose_amb_fec_no_atresia_boxplot.png", plot = glucose_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
glucose_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.), colour = "black") +
  labs(title = "Glucose", x = "Weight Adjusted Fecundity", y = "Glucose (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(glucose_amb_fec_no_atresia_scatterplot)

ggsave(filename = "data-images/glucose_amb_fec_no_atresia_scatterplot.pdf", plot = glucose_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/glucose_amb_fec_no_atresia_scatterplot.pdf", plot = glucose_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image

Glucose Summary Stats

Differences: Between Brood Conditon

# glucose
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

glucose_aov_ambient_brood <- aov(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(glucose_aov_ambient_brood) 
##                           Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition  3  1.218  0.4061   0.837  0.492
## Residuals                 17  8.247  0.4851
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(glucose_aov_ambient_brood) 
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
## $Consensus_Brood_Condition
##                     diff        lwr      upr     p adj
## Normal-Excellent  -0.225 -2.4386100 1.988610 0.9912859
## Low-Excellent     -0.200 -2.4862066 2.086207 0.9943918
## Atresia-Excellent  0.300 -1.7546528 2.354653 0.9751438
## Low-Normal         0.025 -1.4871835 1.537184 0.9999611
## Atresia-Normal     0.525 -0.6070586 1.657059 0.5643929
## Atresia-Low        0.500 -0.7681593 1.768159 0.6821448
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Glu..mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Glu..mmol.L.
## W = 0.85634, p-value = 0.005473
shapiro.test(residuals(aov(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.79895, p-value = 0.0006323
# QQplot:
glucose_ambient_res_qqplot <- ggqqplot(residuals(aov(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Glucose Interactive Residual QQplot",
     subtitle = "residuals(aov(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
                               x = "Glucose Theoretical Quantiles (Predicted values)",
                               y = "Glucose Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(glucose_ambient_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)

# Nonparametric variance test (Levene's test):
leveneTest(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.4505 0.7202
##       17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Glu..mmol.L. by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 3.9713, df = 3, p-value = 0.2646
# Nonparametric Post Hoc: Dunn's test
DunnTest(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                   mean.rank.diff   pval    
## Normal-Excellent      -3.5000000 1.0000    
## Low-Excellent         -2.8333333 1.0000    
## Atresia-Excellent      2.5384615 1.0000    
## Low-Normal             0.6666667 1.0000    
## Atresia-Normal         6.0384615 0.5244    
## Atresia-Low            5.3717949 0.8733    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Glucose Ambient Results

Which brood condition significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartletts test): NA

  • Nonparametric variance (Levene’s test): Met

Differences: Between Fecundity Class

# glucose
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

glucose_aov_ambient_fec <- aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(glucose_aov_ambient_fec)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class  2  1.208  0.6041   1.317  0.293
## Residuals       18  8.257  0.4587
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(glucose_aov_ambient_fec)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                               diff        lwr      upr     p adj
## Low (~1,000s)-High (>50,000) 0.125 -1.0973103 1.347310 0.9632210
## Atresia-High (>50,000)       0.550 -0.4383693 1.538369 0.3518867
## Atresia-Low (~1,000s)        0.425 -0.5633693 1.413369 0.5277964
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Glu..mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Glu..mmol.L.
## W = 0.85634, p-value = 0.005473
shapiro.test(residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.82056, p-value = 0.001376
# QQplot:
glucose_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Glucose Interactive Residual QQplot",
     subtitle = "residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Glucose Theoretical Quantiles (Predicted values)",
                               y = "Glucose Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(glucose_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Glu..mmol.L. by Fecundity_Class
## Bartlett's K-squared = 3.669, df = 2, p-value = 0.1597
# Nonparametric variance test (Levene's test):
leveneTest(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.9391 0.4093
##       18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Glu..mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 4.0755, df = 2, p-value = 0.1303
# Nonparametric Post Hoc: Dunn's test
DunnettTest(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`High (>50,000)`
##                               diff     lwr.ci   upr.ci   pval    
## Low (~1,000s)-High (>50,000) 0.125 -1.0137649 1.263765 0.9470    
## Atresia-High (>50,000)       0.550 -0.3708139 1.470814 0.2758    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Glucose Ambient Results

Which fecundity class significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartlett’s test): Met

Nonparametric variance test (Levene’s test): Met

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# glucose
# Linear Model
glucose_amb_fec_no_atresia_lm <- lm(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(glucose_amb_fec_no_atresia_lm)
## 
## Call:
## lm(formula = Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Fecundity_No_Atresia_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8882 -0.5589 -0.3191  0.2610  1.7458 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                     1.352797   0.563859   2.399
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.002951   0.002957   0.998
##                                                 Pr(>|t|)  
## (Intercept)                                       0.0432 *
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight   0.3475  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8924 on 8 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.1107, Adjusted R-squared:  -0.0004362 
## F-statistic: 0.9961 on 1 and 8 DF,  p-value: 0.3475
glucose_ambient_sim <- simulateResiduals(fittedModel = glucose_amb_fec_no_atresia_lm) 
plot(glucose_ambient_sim)

plotResiduals(glucose_amb_fec_no_atresia_lm)

testDispersion(glucose_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Normality test:
shapiro.test(residuals(glucose_amb_fec_no_atresia_lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(glucose_amb_fec_no_atresia_lm)
## W = 0.86864, p-value = 0.09638
# Transform data: square root 
glucose_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(Glu..mmol.L. = sqrt(Glu..mmol.L.))

# Rerun lm model with square root transformed data
glucose_amb_fec_no_atresia_lm_sqrt <- lm(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt)
summary(glucose_amb_fec_no_atresia_lm_sqrt)
## 
## Call:
## lm(formula = Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = glucose_ambient_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30380 -0.13647 -0.01130  0.07862  0.42723 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.2241905  0.1551624   7.890
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0001917  0.0008838   0.217
##                                                 Pr(>|t|)    
## (Intercept)                                      0.00022 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight  0.83544    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2422 on 6 degrees of freedom
## Multiple R-squared:  0.007783,   Adjusted R-squared:  -0.1576 
## F-statistic: 0.04706 on 1 and 6 DF,  p-value: 0.8354
glucose_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = glucose_amb_fec_no_atresia_lm_sqrt) 
plot(glucose_amb_fec_sqrt_sim)

plotResiduals(glucose_amb_fec_no_atresia_lm_sqrt)

testDispersion(glucose_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
glucose_amb_fec_no_atresia_lm_res <- residuals(glucose_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(glucose_amb_fec_no_atresia_lm_res) 
## 
##  Shapiro-Wilk normality test
## 
## data:  glucose_amb_fec_no_atresia_lm_res
## W = 0.96004, p-value = 0.8104
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt)

Glucose Ambient No Atresia Results

Does a significant correlation exist between weight adj fecundity and parameter?

DHARMA Warning: Data dispersion error

Parametric assumptions: - Normality (Shapiro-Wilk test): Met ish - Variance (Bartlett’s test): NA

Nonparametric variance test (Levene’s test): NA

Residuals Plot (if significant)

# glucose
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

glucose_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Glucose", x = "Weight Adjusted Fecudnity", y = "Glucose (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(glucose_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(filename = "data-images/glucose_ambient_lm_scatterplot.pdf", plot = glucose_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/glucose_ambient_lm_scatterplot.png", plot = glucose_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Na+ : Sodium

Sodium Summary Stats

# sodium
# Summary stats

# Ambient data
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
  summarize(count = n(),
            median = round(median(Na...mmol.L.), 3),
            mean = round(mean(Na...mmol.L.), 3),
            sd = round(sd(Na...mmol.L.), 3),
            cv = round(sd(Na...mmol.L.)/mean(Na...mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
##   Ambient.Or.OAH Consensus_Brood_Condition count median  mean    sd     cv
##   <ord>          <ord>                     <int>  <dbl> <dbl> <dbl>  <dbl>
## 1 Ambient        Excellent                     1   177   177  NA    NA    
## 2 Ambient        Normal                        4   178.  177.  2.87  0.016
## 3 Ambient        Low                           3   175   173.  3.79  0.022
## 4 Ambient        Atresia                      13   178   181. 11.4   0.063
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(Na...mmol.L.), 3),
            mean = round(mean(Na...mmol.L.), 3),
            sd = round(sd(Na...mmol.L.), 3),
            cv = round(sd(Na...mmol.L.)/mean(Na...mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   178.  177.  2.63 0.015
## 2 Ambient        Low (~1,000s)       4   176.  175.  4.19 0.024
## 3 Ambient        Atresia            13   178   181. 11.4  0.063
Amb_Fec_No_Atresia_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(Na...mmol.L.), 3),
            mean = round(mean(Na...mmol.L.), 3),
            sd = round(sd(Na...mmol.L.), 3),
            cv = round(sd(Na...mmol.L.)/mean(Na...mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   178.  177.  2.63 0.015
## 2 Ambient        Low (~1,000s)       4   176.  175.  4.19 0.024

Sodium Boxplot

# Na+ boxplot

# Ambient Samples
sodium_ambient_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = Na...mmol.L.)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Sodium", x = "Parturition Success", y = "Sodium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_ambient_boxplot)

ggsave(filename = "data-images/sodium_ambient_boxplot.pdf", plot = sodium_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/sodium_ambient_boxplot.png", plot = sodium_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
sodium_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_ambient_fecundity_boxplot)

ggsave(filename = "data-images/sodium_ambient_fecundity_boxplot.pdf", plot = sodium_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/sodium_ambient_fecundity_boxplot.png", plot = sodium_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
sodium_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_amb_fec_no_atresia_boxplot)

ggsave(filename = "data-images/sodium_amb_fec_no_atresia_boxplot.pdf", plot = sodium_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/sodium_amb_fec_no_atresia_boxplot.png", plot = sodium_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
sodium_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.), colour = "black") +
  labs(title = "Sodium", x = "Weight Adjusted Fecundity", y = "Sodium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_amb_fec_no_atresia_scatterplot)

ggsave(filename = "data-images/sodium_amb_fec_no_atresia_scatterplot.pdf", plot = sodium_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/sodium_amb_fec_no_atresia_scatterplot.pdf", plot = sodium_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image

Na+ Stat Tests

Differences: Between Brood Conditon

# Na+
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

sodium_aov_ambient_brood <- aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(sodium_aov_ambient_brood) 
##                           Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition  3  176.8   58.95   0.622  0.611
## Residuals                 17 1611.7   94.81
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(sodium_aov_ambient_brood) 
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
## $Consensus_Brood_Condition
##                        diff        lwr      upr     p adj
## Normal-Excellent   0.250000 -30.694633 31.19463 0.9999955
## Low-Excellent     -3.666667 -35.626146 28.29281 0.9875929
## Atresia-Excellent  4.230769 -24.491760 32.95330 0.9745138
## Low-Normal        -3.916667 -25.055875 17.22254 0.9514394
## Atresia-Normal     3.980769 -11.844573 19.80611 0.8897586
## Atresia-Low        7.897436  -9.830494 25.62537 0.5954090
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Na...mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Na...mmol.L.
## W = 0.52676, p-value = 3.625e-07
shapiro.test(residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.54475, p-value = 5.408e-07
# QQplot:
sodium_ambient_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Sodium Residual QQplot",
     subtitle = "residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
                               x = "Sodium Theoretical Quantiles (Predicted values)",
                               y = "Sodium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(sodium_ambient_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)

# Nonparametric variance test (Levene's test):
leveneTest(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   0.194  0.899
##       17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Na...mmol.L. by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 5.5569, df = 3, p-value = 0.1353
# Nonparametric Post Hoc: Dunns
DunnTest(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                   mean.rank.diff   pval    
## Normal-Excellent       3.0000000 1.0000    
## Low-Excellent         -5.5000000 1.0000    
## Atresia-Excellent      3.5769231 1.0000    
## Low-Normal            -8.5000000 0.3481    
## Atresia-Normal         0.5769231 1.0000    
## Atresia-Low            9.0769231 0.1252    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sodium Ambient Results

Which brood condition significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartletts test): NA

  • Nonparametric variance (Levene’s test): Met

  • Nonparametric Post Hoc: Dunn’s Test

  • Atresia-Low: mean.rank.diff = 9.0769231, pval = 0.1252

Differences: Between Fecundity Class

# Na+
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

sodium_aov_ambient_fec <- aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(sodium_aov_ambient_fec)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class  2  156.8   78.38   0.865  0.438
## Residuals       18 1631.8   90.66
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(sodium_aov_ambient_fec)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                   diff       lwr      upr     p adj
## Low (~1,000s)-High (>50,000) -2.000000 -19.18271 15.18271 0.9526465
## Atresia-High (>50,000)        4.480769  -9.41330 18.37484 0.6939179
## Atresia-Low (~1,000s)         6.480769  -7.41330 20.37484 0.4737484
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Na...mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Na...mmol.L.
## W = 0.52676, p-value = 3.625e-07
shapiro.test(residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.56321, p-value = 8.235e-07
# QQplot:
sodium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Sodium Residual QQplot",
     subtitle = "residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Sodium Theoretical Quantiles (Predicted values)",
                               y = "Sodium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(sodium_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Na...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 7.4634, df = 2, p-value = 0.02395
# Nonparametric variance test (Levene's test):
leveneTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.2098 0.8127
##       18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Na...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 2.8608, df = 2, p-value = 0.2392
# Nonparametric Post Hoc: Dunn's test
DunnettTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`High (>50,000)`
##                                   diff     lwr.ci   upr.ci   pval    
## Low (~1,000s)-High (>50,000) -2.000000 -18.008265 14.00826 0.9321    
## Atresia-High (>50,000)        4.480769  -8.463634 17.42517 0.6108    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sodium Ambient Results

Which fecundity class significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartlett’s test): Not Met

Nonparametric variance test (Levene’s test): Met

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# Na+
# Linear Model
sodium_amb_fec_no_atresia_lm <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(sodium_amb_fec_no_atresia_lm)
## 
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Fecundity_No_Atresia_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1021 -0.3010  0.6449  1.9545  3.3374 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     175.27869    2.08086  84.234
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight   0.00558    0.01091   0.511
##                                                 Pr(>|t|)    
## (Intercept)                                      4.4e-13 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.623    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.293 on 8 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.03166,    Adjusted R-squared:  -0.08939 
## F-statistic: 0.2615 on 1 and 8 DF,  p-value: 0.6229
sodium_ambient_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm) 
plot(sodium_ambient_sim)

plotResiduals(sodium_amb_fec_no_atresia_lm)

testDispersion(sodium_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Normailty test
shapiro.test(residuals(sodium_amb_fec_no_atresia_lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(sodium_amb_fec_no_atresia_lm)
## W = 0.86592, p-value = 0.08956
# Transform data: square root 
sodium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(Na...mmol.L. = sqrt(Na...mmol.L.))

# Rerun lm model with square root transformed data
sodium_amb_fec_no_atresia_lm_sqrt <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
summary(sodium_amb_fec_no_atresia_lm_sqrt)
## 
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = sodium_ambient_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25662 -0.03796  0.02115  0.09537  0.12926 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.324e+01  8.923e-02  148.42
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 8.630e-05  5.083e-04    0.17
##                                                 Pr(>|t|)    
## (Intercept)                                     6.31e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.871    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1393 on 6 degrees of freedom
## Multiple R-squared:  0.004782,   Adjusted R-squared:  -0.1611 
## F-statistic: 0.02883 on 1 and 6 DF,  p-value: 0.8708
sodium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm_sqrt) 
plot(sodium_amb_fec_sqrt_sim)
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

plotResiduals(sodium_amb_fec_no_atresia_lm_sqrt)
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

testDispersion(sodium_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
sodium_amb_fec_no_atresia_lm_res <- residuals(sodium_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(sodium_amb_fec_no_atresia_lm_res) 
## 
##  Shapiro-Wilk normality test
## 
## data:  sodium_amb_fec_no_atresia_lm_res
## W = 0.89593, p-value = 0.2654
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)

Sodium Ambient No Atresia Results

Does a significant correlation exist between weight adj fecundity and parameter?

DHARMA Warning: Data dispersion error

Parametric assumptions: - Normality (Shapiro-Wilk test): Not Met - Variance (Bartlett’s test): NA

Nonparametric variance test (Levene’s test): NA

Residuals Plot (if significant)

# Na+
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

sodium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Sodium", x = "Weight Adjusted Fecudnity", y = "Sodium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(filename = "data-images/sodium_ambient_lm_scatterplot.pdf", plot = sodium_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/sodium_ambient_lm_scatterplot.png", plot = sodium_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Cl- : Chloride

Chloride Summary Stats

# Cl-
# Summary stats

# Ambient data
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
  summarize(count = n(),
            median = round(median(Cl...mmol.L.), 3),
            mean = round(mean(Cl...mmol.L.), 3),
            sd = round(sd(Cl...mmol.L.), 3),
            cv = round(sd(Cl...mmol.L.)/mean(Cl...mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
##   Ambient.Or.OAH Consensus_Brood_Condition count median  mean     sd     cv
##   <ord>          <ord>                     <int>  <dbl> <dbl>  <dbl>  <dbl>
## 1 Ambient        Excellent                     1    158  158  NA     NA    
## 2 Ambient        Normal                        4    160  159.  3.86   0.024
## 3 Ambient        Low                           3    156  156.  0.577  0.004
## 4 Ambient        Atresia                      13    158  158.  4.61   0.029
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(Cl...mmol.L.), 3),
            mean = round(mean(Cl...mmol.L.), 3),
            sd = round(sd(Cl...mmol.L.), 3),
            cv = round(sd(Cl...mmol.L.)/mean(Cl...mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   158.  158   2.94 0.019
## 2 Ambient        Low (~1,000s)       4   156   158.  3.70 0.023
## 3 Ambient        Atresia            13   158   158.  4.61 0.029
Amb_Fec_No_Atresia_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(Cl...mmol.L.), 3),
            mean = round(mean(Cl...mmol.L.), 3),
            sd = round(sd(Cl...mmol.L.), 3),
            cv = round(sd(Cl...mmol.L.)/mean(Cl...mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   158.  158   2.94 0.019
## 2 Ambient        Low (~1,000s)       4   156   158.  3.70 0.023

Chloride Boxplot

# Cl- boxplot

# Ambient Samples
chloride_ambient_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = Cl...mmol.L.)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Chloride", x = "Parturition Success", y = "Chloride (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(chloride_ambient_boxplot)

ggsave(filename = "data-images/chloride_ambient_boxplot.pdf", plot = chloride_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/chloride_ambient_boxplot.png", plot = chloride_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
chloride_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Chloride", x = "Fecundity Class", y = "Chloride (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(chloride_ambient_fecundity_boxplot)

ggsave(filename = "data-images/chloride_ambient_fecundity_boxplot.pdf", plot = chloride_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/chloride_ambient_fecundity_boxplot.png", plot = chloride_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
chloride_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Chloride", x = "Fecundity Class", y = "Chloride (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(chloride_amb_fec_no_atresia_boxplot)

ggsave(filename = "data-images/chloride_amb_fec_no_atresia_boxplot.pdf", plot = chloride_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/chloride_amb_fec_no_atresia_boxplot.png", plot = chloride_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
chloride_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.), colour = "black") +
  labs(title = "Chloride", x = "Weight Adjusted Fecundity", y = "Chloride (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(chloride_amb_fec_no_atresia_scatterplot)

ggsave(filename = "data-images/chloride_amb_fec_no_atresia_scatterplot.pdf", plot = chloride_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/chloride_amb_fec_no_atresia_scatterplot.pdf", plot = chloride_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image

Differences: Between Brood Conditon

# Na+
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

sodium_aov_ambient_brood <- aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(sodium_aov_ambient_brood) 
##                           Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition  3  176.8   58.95   0.622  0.611
## Residuals                 17 1611.7   94.81
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(sodium_aov_ambient_brood) 
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
## $Consensus_Brood_Condition
##                        diff        lwr      upr     p adj
## Normal-Excellent   0.250000 -30.694633 31.19463 0.9999955
## Low-Excellent     -3.666667 -35.626146 28.29281 0.9875929
## Atresia-Excellent  4.230769 -24.491760 32.95330 0.9745138
## Low-Normal        -3.916667 -25.055875 17.22254 0.9514394
## Atresia-Normal     3.980769 -11.844573 19.80611 0.8897586
## Atresia-Low        7.897436  -9.830494 25.62537 0.5954090
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Na...mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Na...mmol.L.
## W = 0.52676, p-value = 3.625e-07
# QQplot:
sodium_ambient_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Sodium Interactive Residual QQplot",
     subtitle = "residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
                               x = "Sodium Theoretical Quantiles (Predicted values)",
                               y = "Sodium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(sodium_ambient_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)

# Nonparametric variance test (Levene's test):
leveneTest(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   0.194  0.899
##       17
# Nonparametric stat test: Dunns
DunnTest(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                   mean.rank.diff   pval    
## Normal-Excellent       3.0000000 1.0000    
## Low-Excellent         -5.5000000 1.0000    
## Atresia-Excellent      3.5769231 1.0000    
## Low-Normal            -8.5000000 0.3481    
## Atresia-Normal         0.5769231 1.0000    
## Atresia-Low            9.0769231 0.1252    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sodium Ambient Results

Which brood condition significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartletts test): NA

  • Nonparametric variance (Levene’s test): Met

  • Nonparametric Post Hoc: Dunn’s Test

  • Atresia-Low: mean.rank.diff = 9.0769231, pval = 0.1252

Differences: Between Fecundity Class

# Na+
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

sodium_aov_ambient_fec <- aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(sodium_aov_ambient_fec)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class  2  156.8   78.38   0.865  0.438
## Residuals       18 1631.8   90.66
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(sodium_aov_ambient_fec)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                   diff       lwr      upr     p adj
## Low (~1,000s)-High (>50,000) -2.000000 -19.18271 15.18271 0.9526465
## Atresia-High (>50,000)        4.480769  -9.41330 18.37484 0.6939179
## Atresia-Low (~1,000s)         6.480769  -7.41330 20.37484 0.4737484
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Na...mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Na...mmol.L.
## W = 0.52676, p-value = 3.625e-07
# QQplot:
sodium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Sodium Interactive Residual QQplot",
     subtitle = "residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Sodium Theoretical Quantiles (Predicted values)",
                               y = "Sodium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(sodium_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
#bartlett.test(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

# Nonparametric variance test (Levene's test):
leveneTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.2098 0.8127
##       18
# Nonparametric stat test: Dunns
DunnettTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`High (>50,000)`
##                                   diff     lwr.ci   upr.ci   pval    
## Low (~1,000s)-High (>50,000) -2.000000 -18.008265 14.00826 0.9321    
## Atresia-High (>50,000)        4.480769  -8.463634 17.42517 0.6108    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sodium Ambient Results

Which fecundity class significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartlett’s test): NA

Nonparametric variance test (Levene’s test): Met

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# Na+
# Linear Model
sodium_amb_fec_no_atresia_lm <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(sodium_amb_fec_no_atresia_lm)
## 
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Fecundity_No_Atresia_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1021 -0.3010  0.6449  1.9545  3.3374 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     175.27869    2.08086  84.234
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight   0.00558    0.01091   0.511
##                                                 Pr(>|t|)    
## (Intercept)                                      4.4e-13 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.623    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.293 on 8 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.03166,    Adjusted R-squared:  -0.08939 
## F-statistic: 0.2615 on 1 and 8 DF,  p-value: 0.6229
sodium_ambient_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm) 
plot(sodium_ambient_sim)

plotResiduals(sodium_amb_fec_no_atresia_lm)

testDispersion(sodium_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Transform data: square root 
sodium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(Na...mmol.L. = sqrt(Na...mmol.L.))

# Rerun lm model with square root transformed data
sodium_amb_fec_no_atresia_lm_sqrt <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
summary(sodium_amb_fec_no_atresia_lm_sqrt)
## 
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = sodium_ambient_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25662 -0.03796  0.02115  0.09537  0.12926 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.324e+01  8.923e-02  148.42
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 8.630e-05  5.083e-04    0.17
##                                                 Pr(>|t|)    
## (Intercept)                                     6.31e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.871    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1393 on 6 degrees of freedom
## Multiple R-squared:  0.004782,   Adjusted R-squared:  -0.1611 
## F-statistic: 0.02883 on 1 and 6 DF,  p-value: 0.8708
sodium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm_sqrt) 
plot(sodium_amb_fec_sqrt_sim)
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

plotResiduals(sodium_amb_fec_no_atresia_lm_sqrt)
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

testDispersion(sodium_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
sodium_amb_fec_no_atresia_lm_res <- residuals(sodium_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(sodium_amb_fec_no_atresia_lm_res) 
## 
##  Shapiro-Wilk normality test
## 
## data:  sodium_amb_fec_no_atresia_lm_res
## W = 0.89593, p-value = 0.2654
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)

Sodium Ambient No Atresia Results

Does a significant correlation exist between weight adj fecundity and parameter?

DHARMA Warning: Data dispersion error

Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): NA

Nonparametric variance test (Levene’s test): NA

Residuals Plot (if significant)

# Na+
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

sodium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Sodium", x = "Weight Adjusted Fecudnity", y = "Sodium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(filename = "data-images/sodium_ambient_lm_scatterplot.pdf", plot = sodium_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/sodium_ambient_lm_scatterplot.png", plot = sodium_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Cl- Stat Tests

Differences: Between Brood Condition

# Cl-
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

chloride_aov_ambient_brood <- aov(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(chloride_aov_ambient_brood) 
##                           Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition  3  24.59   8.197   0.463  0.711
## Residuals                 17 300.65  17.685
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(chloride_aov_ambient_brood) 
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
## $Consensus_Brood_Condition
##                         diff        lwr       upr     p adj
## Normal-Excellent   1.2500000 -12.114992 14.614992 0.9931731
## Low-Excellent     -2.3333333 -16.136638 11.469971 0.9623849
## Atresia-Excellent  0.4615385 -11.943726 12.866803 0.9995589
## Low-Normal        -3.5833333 -12.713361  5.546694 0.6851430
## Atresia-Normal    -0.7884615  -7.623430  6.046506 0.9873943
## Atresia-Low        2.7948718  -4.861824 10.451567 0.7302966
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Cl...mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Cl...mmol.L.
## W = 0.96502, p-value = 0.6224
shapiro.test(residuals(aov(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.9502, p-value = 0.3437
# QQplot:
chloride_ambient_res_qqplot <- ggqqplot(residuals(aov(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Chloride Residual QQplot",
     subtitle = "residuals(aov(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
                               x = "Chloride Theoretical Quantiles (Predicted values)",
                               y = "Chloride Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(chloride_ambient_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)

# Nonparametric variance test (Levene's test):
leveneTest(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.3517  0.291
##       17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Cl...mmol.L. by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 2.3847, df = 3, p-value = 0.4965
# Nonparametric Post Hoc: Dunns
DunnTest(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                   mean.rank.diff   pval    
## Normal-Excellent       2.1250000 1.0000    
## Low-Excellent         -4.8333333 1.0000    
## Atresia-Excellent      0.4615385 1.0000    
## Low-Normal            -6.9583333 0.8405    
## Atresia-Normal        -1.6634615 1.0000    
## Atresia-Low            5.2948718 0.9031    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Chloride Ambient Results

Which brood condition significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartletts test): NA

  • Nonparametric variance (Levene’s test): Met

  • Nonparametric Post Hoc: Dunn’s Test

Differences: Between Fecundity Class

# Cl-
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

chloride_aov_ambient_fec <- aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(chloride_aov_ambient_fec)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class  2    3.0   1.504   0.084   0.92
## Residuals       18  322.2  17.902
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(chloride_aov_ambient_fec)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                    diff       lwr      upr     p adj
## Low (~1,000s)-High (>50,000) -0.5000000 -8.135556 7.135556 0.9847331
## Atresia-High (>50,000)        0.4615385 -5.712630 6.635707 0.9801564
## Atresia-Low (~1,000s)         0.9615385 -5.212630 7.135707 0.9170009
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Cl...mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Cl...mmol.L.
## W = 0.96502, p-value = 0.6224
shapiro.test(residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.96555, p-value = 0.6339
# QQplot:
chloride_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Chloride Residual QQplot",
     subtitle = "residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Chloride Theoretical Quantiles (Predicted values)",
                               y = "Chloride Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(chloride_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Cl...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 0.82247, df = 2, p-value = 0.6628
# Nonparametric variance test (Levene's test):
leveneTest(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.5546 0.5838
##       18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Cl...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 0.4005, df = 2, p-value = 0.8185
# Nonparametric Post Hoc: Dunns
DunnettTest(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`High (>50,000)`
##                                    diff    lwr.ci   upr.ci   pval    
## Low (~1,000s)-High (>50,000) -0.5000000 -7.613663 6.613663 0.9777    
## Atresia-High (>50,000)        0.4615385 -5.290623 6.213700 0.9711    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sodium Ambient Results

Which fecundity class significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): Met

Nonparametric variance test (Levene’s test): Met

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# Cl-
# Linear Model
chloride_amb_fec_no_atresia_lm <- lm(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(chloride_amb_fec_no_atresia_lm)
## 
## Call:
## lm(formula = Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Fecundity_No_Atresia_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1972 -2.1453 -0.3044  2.3202  4.8677 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.581e+02  2.042e+00  77.415
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 7.028e-04  1.071e-02   0.066
##                                                 Pr(>|t|)    
## (Intercept)                                     8.64e-13 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.949    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.232 on 8 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.0005381,  Adjusted R-squared:  -0.1244 
## F-statistic: 0.004307 on 1 and 8 DF,  p-value: 0.9493
chloride_ambient_sim <- simulateResiduals(fittedModel = chloride_amb_fec_no_atresia_lm) 
plot(chloride_ambient_sim)

plotResiduals(chloride_amb_fec_no_atresia_lm)

testDispersion(chloride_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Shapiro Test
shapiro.test(residuals(chloride_amb_fec_no_atresia_lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(chloride_amb_fec_no_atresia_lm)
## W = 0.95118, p-value = 0.6824
# Transform data: square root 
chloride_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(Cl...mmol.L. = sqrt(Cl...mmol.L.))

# Rerun lm model with square root transformed data
chloride_amb_fec_no_atresia_lm_sqrt <- lm(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt)
summary(chloride_amb_fec_no_atresia_lm_sqrt)
## 
## Call:
## lm(formula = Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = chloride_ambient_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14856 -0.08721 -0.02279  0.07258  0.20200 
## 
## Coefficients:
##                                                   Estimate Std. Error t value
## (Intercept)                                      1.257e+01  8.518e-02 147.566
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -7.493e-05  4.852e-04  -0.154
##                                                 Pr(>|t|)    
## (Intercept)                                     6.53e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.882    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.133 on 6 degrees of freedom
## Multiple R-squared:  0.003959,   Adjusted R-squared:  -0.162 
## F-statistic: 0.02385 on 1 and 6 DF,  p-value: 0.8823
chloride_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = chloride_amb_fec_no_atresia_lm_sqrt) 
plot(chloride_amb_fec_sqrt_sim)
## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.

plotResiduals(chloride_amb_fec_no_atresia_lm_sqrt)
## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.

testDispersion(chloride_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
chloride_amb_fec_no_atresia_lm_res <- residuals(chloride_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(chloride_amb_fec_no_atresia_lm_res) 
## 
##  Shapiro-Wilk normality test
## 
## data:  chloride_amb_fec_no_atresia_lm_res
## W = 0.94043, p-value = 0.6153
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt)

Chloride Ambient No Atresia Results

Does a significant correlation exist between weight adj fecundity and parameter?

DHARMA Warning: Data dispersion error

Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): NA

Nonparametric variance test (Levene’s test): NA

Residuals Plot (if significant)

# Cl-
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

chloride_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Cl...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Chloride", x = "Weight Adjusted Fecudnity", y = "Chloride (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(chloride_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(filename = "data-images/chloride_ambient_lm_scatterplot.pdf", plot = chloride_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/chloride_ambient_lm_scatterplot.png", plot = chloride_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

K+ : Potassium

Potassium Summary Stats

# K+
# Summary stats

# Ambient data
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
  summarize(count = n(),
            median = round(median(K...mmol.L.), 3),
            mean = round(mean(K...mmol.L.), 3),
            sd = round(sd(K...mmol.L.), 3),
            cv = round(sd(K...mmol.L.)/mean(K...mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
##   Ambient.Or.OAH Consensus_Brood_Condition count median  mean     sd     cv
##   <ord>          <ord>                     <int>  <dbl> <dbl>  <dbl>  <dbl>
## 1 Ambient        Excellent                     1    2.7  2.7  NA     NA    
## 2 Ambient        Normal                        4    2.7  2.7   0.082  0.03 
## 3 Ambient        Low                           3    2.7  2.73  0.252  0.092
## 4 Ambient        Atresia                      13    2.7  2.86  0.304  0.106
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(K...mmol.L.), 3),
            mean = round(mean(K...mmol.L.), 3),
            sd = round(sd(K...mmol.L.), 3),
            cv = round(sd(K...mmol.L.)/mean(K...mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4    2.7  2.7  0.082 0.03 
## 2 Ambient        Low (~1,000s)       4    2.7  2.72 0.206 0.076
## 3 Ambient        Atresia            13    2.7  2.86 0.304 0.106
Amb_Fec_No_Atresia_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(K...mmol.L.), 3),
            mean = round(mean(K...mmol.L.), 3),
            sd = round(sd(K...mmol.L.), 3),
            cv = round(sd(K...mmol.L.)/mean(K...mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4    2.7  2.7  0.082 0.03 
## 2 Ambient        Low (~1,000s)       4    2.7  2.72 0.206 0.076

Potassium Boxplot

# K+ boxplot

# Ambient Samples
potassium_ambient_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = K...mmol.L.)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = K...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Potassium", x = "Parturition Success", y = "Potassium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(potassium_ambient_boxplot)

ggsave(filename = "data-images/potassium_ambient_boxplot.pdf", plot = potassium_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/potassium_ambient_boxplot.png", plot = potassium_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
potassium_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = K...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = K...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Potassium", x = "Fecundity Class", y = "Potassium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(potassium_ambient_fecundity_boxplot)

ggsave(filename = "data-images/potassium_ambient_fecundity_boxplot.pdf", plot = potassium_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/potassium_ambient_fecundity_boxplot.png", plot = potassium_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
potassium_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = K...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = K...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Potassium", x = "Fecundity Class", y = "Potassium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(potassium_amb_fec_no_atresia_boxplot)

ggsave(filename = "data-images/potassium_amb_fec_no_atresia_boxplot.pdf", plot = potassium_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/potassium_amb_fec_no_atresia_boxplot.png", plot = potassium_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
potassium_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.), colour = "black") +
  labs(title = "Potassium", x = "Weight Adjusted Fecundity", y = "Potassium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(potassium_amb_fec_no_atresia_scatterplot)

ggsave(filename = "data-images/potassium_amb_fec_no_atresia_scatterplot.pdf", plot = potassium_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/potassium_amb_fec_no_atresia_scatterplot.pdf", plot = potassium_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image

K+ Stat Tests

Differences: Between Brood Condition

# K+
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

potassium_aov_ambient_brood <- aov(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(potassium_aov_ambient_brood) 
##                           Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition  3 0.1121 0.03736   0.505  0.684
## Residuals                 17 1.2574 0.07397
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(potassium_aov_ambient_brood) 
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
## $Consensus_Brood_Condition
##                         diff        lwr       upr     p adj
## Normal-Excellent  0.00000000 -0.8643366 0.8643366 1.0000000
## Low-Excellent     0.03333333 -0.8593496 0.9260163 0.9995541
## Atresia-Excellent 0.16153846 -0.6407309 0.9638078 0.9389808
## Low-Normal        0.03333333 -0.5571209 0.6237876 0.9984688
## Atresia-Normal    0.16153846 -0.2804904 0.6035674 0.7296107
## Atresia-Low       0.12820513 -0.3669663 0.6233765 0.8812696
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$K...mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$K...mmol.L.
## W = 0.81347, p-value = 0.001062
shapiro.test(residuals(aov(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.886, p-value = 0.0189
# QQplot:
potassium_ambient_res_qqplot <- ggqqplot(residuals(aov(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Potassium Residual QQplot",
     subtitle = "residuals(aov(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
                               x = "Potassium Theoretical Quantiles (Predicted values)",
                               y = "Potassium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(potassium_ambient_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)

# Nonparametric variance test (Levene's test):
leveneTest(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.6208 0.6111
##       17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  K...mmol.L. by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 1.2134, df = 3, p-value = 0.7498
# Nonparametric Post Hoc: Dunn's test
DunnTest(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                   mean.rank.diff   pval    
## Normal-Excellent        0.000000 1.0000    
## Low-Excellent           0.500000 1.0000    
## Atresia-Excellent       3.115385 1.0000    
## Low-Normal              0.500000 1.0000    
## Atresia-Normal          3.115385 1.0000    
## Atresia-Low             2.615385 1.0000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Potassium Ambient Results

Which brood condition significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Not Met - Variance (Bartletts test): NA

  • Nonparametric variance (Levene’s test): Met

  • Nonparametric Post Hoc: Dunn’s Test

Differences: Between Fecundity Class

# K+
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

potassium_aov_ambient_fec <- aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(potassium_aov_ambient_fec)
##                 Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class  2 0.1113 0.05563   0.796  0.466
## Residuals       18 1.2583 0.06990
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(potassium_aov_ambient_fec)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                   diff        lwr       upr     p adj
## Low (~1,000s)-High (>50,000) 0.0250000 -0.4521380 0.5021380 0.9901955
## Atresia-High (>50,000)       0.1615385 -0.2242789 0.5473558 0.5449631
## Atresia-Low (~1,000s)        0.1365385 -0.2492789 0.5223558 0.6452587
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$K...mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$K...mmol.L.
## W = 0.81347, p-value = 0.001062
shapiro.test(residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.88368, p-value = 0.0171
# QQplot:
potassium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Potassium Residual QQplot",
     subtitle = "residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Potassium Theoretical Quantiles (Predicted values)",
                               y = "Potassium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(potassium_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  K...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 4.6368, df = 2, p-value = 0.09843
# Nonparametric variance test (Levene's test):
leveneTest(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.7708 0.4773
##       18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  K...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 1.2081, df = 2, p-value = 0.5466
# Nonparametric stat test: Kruskall Wallis
kruskal.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  K...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 1.2081, df = 2, p-value = 0.5466
# Nonparametric Post Hoc: Dunns test
DunnettTest(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`High (>50,000)`
##                                   diff     lwr.ci    upr.ci   pval    
## Low (~1,000s)-High (>50,000) 0.0250000 -0.4195254 0.4695254 0.9857    
## Atresia-High (>50,000)       0.1615385 -0.1979081 0.5209850 0.4550    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Potassium Ambient Results

Which fecundity class significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Not Met - Variance (Bartlett’s test): Kinda Met

Nonparametric variance test (Levene’s test): Met

Nonparametric Post Hoc: Kruskall Wallis test & Dunns test

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# K+
# Linear Model
potassium_amb_fec_no_atresia_lm <- lm(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(potassium_amb_fec_no_atresia_lm)
## 
## Call:
## lm(formula = K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Fecundity_No_Atresia_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.21166 -0.11568 -0.04065  0.00147  0.39801 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     2.6473078  0.1245207  21.260
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0007432  0.0006530   1.138
##                                                 Pr(>|t|)    
## (Intercept)                                     2.52e-08 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.288    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1971 on 8 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.1393, Adjusted R-squared:  0.03177 
## F-statistic: 1.295 on 1 and 8 DF,  p-value: 0.288
potassium_ambient_sim <- simulateResiduals(fittedModel = potassium_amb_fec_no_atresia_lm) 
plot(potassium_ambient_sim)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully

plotResiduals(potassium_amb_fec_no_atresia_lm)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully

testDispersion(potassium_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Shapiro Test
shapiro.test(residuals(potassium_amb_fec_no_atresia_lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(potassium_amb_fec_no_atresia_lm)
## W = 0.86151, p-value = 0.0795
# Transform data: square root 
potassium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(K...mmol.L. = sqrt(K...mmol.L.))

# Rerun lm model with square root transformed data
potassium_amb_fec_no_atresia_lm_sqrt <- lm(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)
summary(potassium_amb_fec_no_atresia_lm_sqrt)
## 
## Call:
## lm(formula = K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = potassium_ambient_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.043790 -0.030608 -0.001939  0.010984  0.085420 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.6239747  0.0283420  57.299
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0001535  0.0001614   0.951
##                                                 Pr(>|t|)    
## (Intercept)                                      1.9e-09 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.378    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04424 on 6 degrees of freedom
## Multiple R-squared:  0.131,  Adjusted R-squared:  -0.01382 
## F-statistic: 0.9046 on 1 and 6 DF,  p-value: 0.3783
potassium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = potassium_amb_fec_no_atresia_lm_sqrt) 
plot(potassium_amb_fec_sqrt_sim)
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## qu = 0.75, log(sigma) = -3.903504 : outer Newton did not converge fully.

plotResiduals(potassium_amb_fec_no_atresia_lm_sqrt)
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## qu = 0.75, log(sigma) = -3.903504 : outer Newton did not converge fully.

testDispersion(potassium_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
potassium_amb_fec_no_atresia_lm_res <- residuals(potassium_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(potassium_amb_fec_no_atresia_lm_res) 
## 
##  Shapiro-Wilk normality test
## 
## data:  potassium_amb_fec_no_atresia_lm_res
## W = 0.88746, p-value = 0.2216
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
# bartlett.test(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)

Potassium Ambient No Atresia Results

Does a significant correlation exist between weight adj fecundity and parameter?

DHARMA Warning: Data dispersion error

Parametric assumptions: - Normality (Shapiro-Wilk test): Kinda Met - Variance (Bartlett’s test): NA

Nonparametric variance test (Levene’s test): NA

Residuals Plot (if significant)

# K+
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

potassium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Potassium", x = "Weight Adjusted Fecudnity", y = "Potassium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(potassium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(filename = "data-images/potassium_ambient_lm_scatterplot.pdf", plot = potassium_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/potassium_ambient_lm_scatterplot.png", plot = potassium_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Ca+2 : Calcium

Calcium Summary Stats

# Ca+2
# Summary stats

# Ambient data
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
  summarize(count = n(),
            median = round(median(Ca....mmol.L.), 3),
            mean = round(mean(Ca....mmol.L.), 3),
            sd = round(sd(Ca....mmol.L.), 3),
            cv = round(sd(Ca....mmol.L.)/mean(Ca....mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
##   Ambient.Or.OAH Consensus_Brood_Condition count median  mean     sd     cv
##   <ord>          <ord>                     <int>  <dbl> <dbl>  <dbl>  <dbl>
## 1 Ambient        Excellent                     1   1.47  1.47 NA     NA    
## 2 Ambient        Normal                        4   1.42  1.43  0.05   0.035
## 3 Ambient        Low                           3   1.43  1.42  0.036  0.025
## 4 Ambient        Atresia                      13   1.24  1.27  0.114  0.09
Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(Ca....mmol.L.), 3),
            mean = round(mean(Ca....mmol.L.), 3),
            sd = round(sd(Ca....mmol.L.), 3),
            cv = round(sd(Ca....mmol.L.)/mean(Ca....mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   1.45  1.45 0.044 0.03 
## 2 Ambient        Low (~1,000s)       4   1.41  1.41 0.033 0.023
## 3 Ambient        Atresia            13   1.24  1.27 0.114 0.09
Amb_Fec_No_Atresia_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarize(count = n(),
            median = round(median(Ca....mmol.L.), 3),
            mean = round(mean(Ca....mmol.L.), 3),
            sd = round(sd(Ca....mmol.L.), 3),
            cv = round(sd(Ca....mmol.L.)/mean(Ca....mmol.L.), 3)) %>%
  ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
##   Ambient.Or.OAH Fecundity_Class count median  mean    sd    cv
##   <ord>          <ord>           <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Ambient        High (>50,000)      4   1.45  1.45 0.044 0.03 
## 2 Ambient        Low (~1,000s)       4   1.41  1.41 0.033 0.023

Calcium Boxplot

# Ca+2 boxplot

# Ambient Samples
calcium_ambient_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = Ca....mmol.L.)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Calcium", x = "Parturition Success", y = "Calcium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(calcium_ambient_boxplot)

ggsave(filename = "data-images/calcium_ambient_boxplot.pdf", plot = calcium_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/calcium_ambient_boxplot.png", plot = calcium_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
calcium_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Ca....mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Calcium", x = "Fecundity Class", y = "Calcium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(calcium_ambient_fecundity_boxplot)

ggsave(filename = "data-images/calcium_ambient_fecundity_boxplot.pdf", plot = calcium_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/calcium_ambient_fecundity_boxplot.png", plot = calcium_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
calcium_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Ca....mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Calcium", x = "Fecundity Class", y = "Calcium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(calcium_amb_fec_no_atresia_boxplot)

ggsave(filename = "data-images/calcium_amb_fec_no_atresia_boxplot.pdf", plot = calcium_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/calcium_amb_fec_no_atresia_boxplot.png", plot = calcium_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
calcium_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.), colour = "black") +
  labs(title = "Calcium", x = "Weight Adjusted Fecundity", y = "Calcium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(calcium_amb_fec_no_atresia_scatterplot)

ggsave(filename = "data-images/calcium_amb_fec_no_atresia_scatterplot.pdf", plot = calcium_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/calcium_amb_fec_no_atresia_scatterplot.pdf", plot = calcium_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image

Ca+2 Stat Tests

Differences: Between Brood Condition

# Ca+2
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

calcium_aov_ambient_brood <- aov(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(calcium_aov_ambient_brood) 
##                           Df Sum Sq Mean Sq F value Pr(>F)  
## Consensus_Brood_Condition  3 0.1319 0.04396   4.527 0.0165 *
## Residuals                 17 0.1651 0.00971                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(calcium_aov_ambient_brood) 
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
## $Consensus_Brood_Condition
##                         diff        lwr           upr     p adj
## Normal-Excellent  -0.0400000 -0.3531867  0.2731867446 0.9830709
## Low-Excellent     -0.0500000 -0.3734579  0.2734578789 0.9707679
## Atresia-Excellent -0.2007692 -0.4914663  0.0899278749 0.2401278
## Low-Normal        -0.0100000 -0.2239473  0.2039472768 0.9991279
## Atresia-Normal    -0.1607692 -0.3209355 -0.0006029264 0.0489650
## Atresia-Low       -0.1507692 -0.3301914  0.0286529182 0.1173091
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Ca....mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Ca....mmol.L.
## W = 0.9463, p-value = 0.2894
shapiro.test(residuals(aov(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.91568, p-value = 0.07117
# QQplot:
calcium_ambient_res_qqplot <- ggqqplot(residuals(aov(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Calcium Residual QQplot",
     subtitle = "residuals(aov(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
                               x = "Calcium Theoretical Quantiles (Predicted values)",
                               y = "Calcium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(calcium_ambient_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)

# Nonparametric variance test (Levene's test):
leveneTest(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.9908 0.4207
##       17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Ca....mmol.L. by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 9.0798, df = 3, p-value = 0.02825
# Nonparametric Post Hoc: Dunn's test
DunnTest(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                   mean.rank.diff   pval    
## Normal-Excellent       -3.125000 1.0000    
## Low-Excellent          -3.500000 1.0000    
## Atresia-Excellent     -11.153846 0.3321    
## Low-Normal             -0.375000 1.0000    
## Atresia-Normal         -8.028846 0.1412    
## Atresia-Low            -7.653846 0.2698    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Potassium Ambient Results

Which brood condition significantly different (if any)? - Atresia-Normal : (Tukey HSD: p adj = 0.0489650) - Atresia-Low : (Tueky HSD: p adj = 0.1173091)

Parametric assumptions: - Normality (Shapiro-Wilk test): Maybe Met - Variance (Bartletts test): NA

Differences: Between Fecundity Class

# Ca+2
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

calcium_aov_ambient_fec <- aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(calcium_aov_ambient_fec)
##                 Df Sum Sq Mean Sq F value  Pr(>F)   
## Fecundity_Class  2 0.1328 0.06641   7.281 0.00482 **
## Residuals       18 0.1642 0.00912                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(potassium_aov_ambient_fec)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                   diff        lwr       upr     p adj
## Low (~1,000s)-High (>50,000) 0.0250000 -0.4521380 0.5021380 0.9901955
## Atresia-High (>50,000)       0.1615385 -0.2242789 0.5473558 0.5449631
## Atresia-Low (~1,000s)        0.1365385 -0.2492789 0.5223558 0.6452587
# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Ca....mmol.L.)
## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Ca....mmol.L.
## W = 0.9463, p-value = 0.2894
shapiro.test(residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.91863, p-value = 0.08147
# QQplot:
calcium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Calcium Residual QQplot",
     subtitle = "residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Calcium Theoretical Quantiles (Predicted values)",
                               y = "Calcium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(calcium_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ca....mmol.L. by Fecundity_Class
## Bartlett's K-squared = 6.1296, df = 2, p-value = 0.04666
# Nonparametric variance test (Levene's test):
leveneTest(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.2704 0.3047
##       18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Ca....mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 9.153, df = 2, p-value = 0.01029
# Nonparametric Post Hoc: Dunns test
DunnettTest(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`High (>50,000)`
##                                    diff     lwr.ci     upr.ci   pval    
## Low (~1,000s)-High (>50,000) -0.0375000 -0.1980658  0.1230658 0.7893    
## Atresia-High (>50,000)       -0.1807692 -0.3106040 -0.0509345 0.0071 ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Potassium Ambient Results

Which fecundity class significantly different (if any)?

Parametric assumptions: - Normality (Shapiro-Wilk test): Maybe Met - Variance (Bartlett’s test): Not Met

Nonparametric variance test (Levene’s test): Met

Nonparametric Post Hoc: Kruskall Wallis test & Dunns test - Atresia-High (>50,000) : (Dunn’s test: pval = 0.0071 **)

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# Ca+2
# Linear Model
calcium_amb_fec_no_atresia_lm <- lm(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(calcium_amb_fec_no_atresia_lm)
## 
## Call:
## lm(formula = Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Fecundity_No_Atresia_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.19339 -0.01787  0.01808  0.04987  0.10485 
## 
## Coefficients:
##                                                   Estimate Std. Error t value
## (Intercept)                                      1.401e+00  5.882e-02  23.823
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -3.752e-05  3.084e-04  -0.122
##                                                 Pr(>|t|)    
## (Intercept)                                     1.03e-08 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.906    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09309 on 8 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.001846,   Adjusted R-squared:  -0.1229 
## F-statistic: 0.01479 on 1 and 8 DF,  p-value: 0.9062
calcium_ambient_sim <- simulateResiduals(fittedModel = calcium_amb_fec_no_atresia_lm) 
plot(calcium_ambient_sim)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully

plotResiduals(calcium_amb_fec_no_atresia_lm)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully

testDispersion(calcium_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Shapiro Test
shapiro.test(residuals(calcium_amb_fec_no_atresia_lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(calcium_amb_fec_no_atresia_lm)
## W = 0.9084, p-value = 0.2702
# Transform data: square root 
calcium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(Ca....mmol.L. = sqrt(Ca....mmol.L.))

# Rerun lm model with square root transformed data
calcium_amb_fec_no_atresia_lm_sqrt <- lm(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)
summary(calcium_amb_fec_no_atresia_lm_sqrt)
## 
## Call:
## lm(formula = Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = calcium_ambient_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.009302 -0.007894 -0.006232  0.004396  0.026830 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.180e+00  8.621e-03 136.816
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 1.143e-04  4.911e-05   2.328
##                                                 Pr(>|t|)    
## (Intercept)                                     1.03e-11 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight   0.0588 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01346 on 6 degrees of freedom
## Multiple R-squared:  0.4745, Adjusted R-squared:  0.3869 
## F-statistic: 5.418 on 1 and 6 DF,  p-value: 0.05882
calcium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = calcium_amb_fec_no_atresia_lm_sqrt) 
plot(calcium_amb_fec_sqrt_sim)
## qu = 0.25, log(sigma) = -3.52446 : outer Newton did not converge fully.

plotResiduals(calcium_amb_fec_no_atresia_lm_sqrt)
## qu = 0.25, log(sigma) = -3.52446 : outer Newton did not converge fully.

testDispersion(calcium_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
calcium_amb_fec_no_atresia_lm_res <- residuals(calcium_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(calcium_amb_fec_no_atresia_lm_res) 
## 
##  Shapiro-Wilk normality test
## 
## data:  calcium_amb_fec_no_atresia_lm_res
## W = 0.77335, p-value = 0.01475
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
# bartlett.test(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)

Potassium Ambient No Atresia Results

Does a significant correlation exist between weight adj fecundity and parameter?

Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): NA

Nonparametric variance test (Levene’s test): NA

Residuals Plot (if significant)

# Ca+2
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

calcium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Calcium", x = "Weight Adjusted Fecudnity", y = "Calcium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(calcium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(filename = "data-images/calcium_ambient_lm_scatterplot.pdf", plot = calcium_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/calcium_ambient_lm_scatterplot.png", plot = calcium_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Combined plots

# Brood Condition Boxplots

ph_amb_bc <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = pH)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = pH), alpha = 0.5, colour = "black") +
  labs(title = "pH", x = "Parturition Success", y = "pH") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

hct_amb_bc <-  ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = Hct....)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = Hct....), alpha = 0.5, colour = "black") +
  labs(title = "Hematocrit", x = "Parturition Success", y = "Hematocrit (%)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
        
        
glu_amb_bc <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = Glu..mmol.L.)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Glucose", x = "Parturition Success", y = "Glucose (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

na_amb_bc <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

cl_amb_bc <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = Cl...mmol.L.)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Chloride", x = "Parturition Success", y = "Chloride (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

k_amb_bc <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = K...mmol.L.)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = K...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Potassium", x = "Parturition Success", y = "Potassium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

ca_amb_bc <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Consensus_Brood_Condition, y = Ca....mmol.L.)) +
  geom_point(aes(x = Consensus_Brood_Condition, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Calcium", x = "Parturition Success", y = "Calcium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

brood_condition_boxplots <- ph_amb_bc + hct_amb_bc + glu_amb_bc + na_amb_bc + cl_amb_bc + k_amb_bc + ca_amb_bc + plot_layout(ncol = 3, guides = "collect") 

brood_condition_boxplots

# Fecundity class boxplots 

ph_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = pH)) +
  geom_point(aes(x = Fecundity_Class, y = pH), alpha = 0.5, colour = "black") +
  labs(title = "pH", x = "Fecundity Class", y = "pH") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

hct_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Hct....)) +
  geom_point(aes(x = Fecundity_Class, y = Hct....), alpha = 0.5, colour = "black") +
  labs(title = "Hematocrit", x = "Fecundity Class", y = "Hematocrit (%)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  

glu_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Glu..mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Glucose", x = "Fecundity Class", y = "Glucose (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  
  
na_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  

cl_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Chloride", x = "Fecundity Class", y = "Chloride (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  

k_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = K...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = K...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Potassium", x = "Fecundity Class", y = "Potassium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  
  
ca_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Ca....mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Calcium", x = "Fecundity Class", y = "Calcium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  

fec_class_boxplots <- ph_fec_class + hct_fec_class + glu_fec_class + na_fec_class + cl_fec_class + k_fec_class + ca_fec_class + plot_layout(guides = "collect")

fec_class_boxplots

All lm plots

# lm plots
ph_am_lm <- ggplot(data = pH_ambient_sqrt, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Weight Adjusted Fecudnity", y = "Sqrt pH", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

hct_am_lm <- ggplot(data = hct_ambient_sqrt, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Hematocrit", x = "Weight Adjusted Fecudnity", y = "Sqrt Hematocrit (%)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

glucose_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Glucose", x = "Weight Adjusted Fecudnity", y = "Glucose (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

sodium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Sodium", x = "Weight Adjusted Fecudnity", y = "Sodium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

chloride_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Cl...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Chloride", x = "Weight Adjusted Fecudnity", y = "Chloride (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

potassium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Potassium", x = "Weight Adjusted Fecudnity", y = "Potassium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

calcium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Calcium", x = "Weight Adjusted Fecudnity", y = "Calcium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
        
ambient_lm_plots <- ph_am_lm + hct_am_lm + glucose_am_lm + sodium_am_lm + chloride_am_lm + potassium_am_lm + calcium_am_lm + plot_layout(guides = "collect")

ambient_lm_plots
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).